Spectral properties of strongly correlated bosons in two-dimensional optical lattices 
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Spectral properties of the two-dimensional Bose-Hubbard model, which emulates ultracold gases of 
atoms confined in optical lattices, are investigated by means of the variational cluster approach. The 
phase boundary of the quantum phase transition from Mott phase to superfluid phase is calculated 
and compared to recent work. Moreover the single-particle spectral functions in both the first 
and the second Mott lobe are presented and the corresponding densities of states and momentum 
distributions are evaluated. A qualitatively similar intensity distribution of the spectral weight can 
be observed for spectral functions in the first and the second Mott lobe. 
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I. INTRODUCTION 

Pioneering experiments on ultracold gases of atoms 
trapped in optical lattices allowed for a direct observation 
of quantum many-body phenomena, such as the quantum 
phase transition from Mott phase to superfluid phase 
Optical lattices are realized by countcrpropagating laser 
beams, which form a periodic potential^ The bosonic 
particles located on the optical lattice gain kinetic en- 
ergy when tunneling through the potential wells of neigh- 
boring sites of the periodic potential and they exhibit a 
repulsive interaction when a lattice site is occupied by 
more than one atom. A condensate of ultracold atoms 
can be driven from superfluid phase to Mott phase by 
gradually increasing the intensity of the laser beams. The 
potential wells of the optical lattice are shallow for low 
laser-beam intensity. Thus the bosonic particles can over- 
come the barrier easily and are delocalized on the whole 
lattice. However, for large intensity of the laser beams 
the potential wells are deep and there is little proba- 
bility for the atoms to tunnel from one lattice site to 
another. This physical behavior can be described by the 
Bose-Hubbard (BH) model^ provided the gas of ultracold 
atoms is cooled such that only the lowest Bloch band of 
the periodic potential has to be taken into account^ The 
ground state of the BH model is superfluid when the local 
on-site repulsion between the atoms is small in compari- 
son to the nearest-neighbor hopping strength whereas it 
is a Mott state for integer particle density and large on- 
site repulsion compared to the hopping strength. Due 
to these characteristics of the BH model the depth of 
the potential wells in optical lattices can be associated 
directly with the ratio of the on-site repulsion and the 
hopping strength. Ultracold atoms confined in optical 
lattices provide a very clean experimental realization of 
a strongly correlated many-body problem and the inter- 
nal physical processes are well understood in comparison 
to conventional condensed-matter systems. There is large 
experimental control over the system parameters, such as 
the particle number, lattice size, and depth of the poten- 
tial wells. Furthermore the sites of the optical lattice can 
be addressed individually due to the mcsoscopic scale of 



the latticed 

The quantum phase transition from Mott phase to su- 
perfluid phase has been first observed experimentally for 
ultracold rubidium atoms trapped in a three-dimensional 
optical lattice^ and subsequently as well in optical lat- 
tices of two dimensions The corresponding theoretical 
model, the two-dimensional (2D) BH model, has already 
been investigated to some detail in literature. The phase 
diagram, which describes the quantum phase transition 
from Mott phase to superfluid phase, has been investi- 
gated thoroughly at the mean-field level (possibly includ- 
ing Gaussian- fluctuation corrections) More accu- 
rate results for the phase diagram from quantum Monte 
Carloi^ (QMC) simulations, variational approaches ; 15 : 16 
and strong-coupling perturbation theory^ 7 - - — are also 
available. The phase diagram for arbitrary integer fill- 
ings has been obtained recently using the so-called di- 
agrammatic process chain approacln 21 ' 22 Spectral func- 
tions of the two-dimensional BH model have been evalu- 
ated within a strong-coupling approac h 19 ' 23 and a varia- 
tional mean field approach^ 

In the present paper we evaluate the border of the 
quantum phase transition from Mott phase to superfluid 
phase for the first two Mott lobes by means of the vari- 
ational cluster approach (VCA)^ and show that this 
method provides quite accurately the boundaries of the 
Mott phase, as compared with more demanding QMC 
simulations and perturbative expansions. In addition, 
we study in detail the spectral functions of the two- 
dimensional BH model in both the first and the second 
Mott lobe, which require computing the Green's function 
in real frequency domain. We also present the densities 
of states and momentum distributions corresponding to 
the spectral functions. Finally, as a technical point, we 
present an extension of the so-called Q-matrix formalism, 
which has been originally proposed for fcrmionic (anti- 
commutator) Green's functions ) 26 ' 27 to bosonic (commu- 
tator) Green's functions.— As we show below, this ex- 
tension is nontrivial due to the nonunitary nature of the 
Bogoliubov transformation for bosonic particles. 

This paper is organized as follows. In Sec. |H] the BH 
model is introduced. Section Mil contains a short review 
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on the variational cluster approach and the extension of 
the Q-matrix formalism. Section IIVI is devoted to the 
spectral properties of the BH model in two dimensions. 
Here the phase diagram, spectral functions, densities of 
states, and momentum distributions are presented. Fi- 
nally, we summarize and conclude our findings in Sec.fVl 

II. MODEL 

The (grand-canonical) Hamiltonian of the BH model^ 
is given by 

H=- t Y. ( 6 1 b J + H ' c -) + 2 ("< - 1) - 

(1) 

where t is the nearest-neighbor hopping strength, U is the 
local on-site repulsion, and [l is the chemical potential. 
The angle brackets in the first part of the Hamiltonian 
specify to sum over pairs of nearest neighbors (each pair 
counted once). The operator b\ creates a particle at lat- 
tice site i whereas b% annihilates a particle at site i. The 
total particle number 

i i 

is conserved, since [H, N p ] = 0. The particles of the BH 
model are of bosonic character and thus the commutation 
relation [bi, ftt] = 5ij is satisfied. The first term of the 
Hamiltonian models the hopping of a particle from lat- 
tice site j to lattice site i. The second part describes the 
local on-site repulsion, which remains zero when a lat- 
tice site is unoccupied or occupied by only one particle. 
However, it increases proportional to U for each addition- 
ally added particle. We consider the on-site repulsion U 
as unit of energy. The third part of the Hamiltonian is 
necessary to perform calculations in grand-canonical en- 
semble, where the chemical potential fj, controls the total 
particle number of the system. 

III. METHOD 
A. Variational cluster approach 

We use VCA3£ to evaluate the phase diagram and the 
spectral functions of the 2D BH model. VCA is a vari- 
ational extension of the cluster perturbation thcor y 29 i 30 
and is based on the self-energy functional approach (SFA) 
which has been originally proposed for fcrmionic systems 
by M. Potthoffi2ii22 VCA has been extended to bosonic 
systems as welli^ 3 . 

The SFA is based on the fact that Dyson's equation 
for the exact Green's function is recovered at the sta- 
tionary point of the grand potential f2[Z] considered as 
a functional of the self-energy Z. Thus Z corresponds, 
at the stationary point, to the real physical self-energy 



The self-energy functional fi[Z] cannot be evaluated di- 
rectly as it contains the Legendre transform F[Y] of the 
Luttinger-Ward functional ^ 1 ' 34 However, the functional 
F[L] just depends on the interaction term of the Hamil- 
tonian, i.e., on the second term of Eq. (|Tj), and is thus 
equivalent for all Hamiltonians which share a common 
interaction part. Due to this property F[T] can be elim- 
inated from the expression of the self-energy functional 
f2[Z]. For this purpose an exactly solvable, so-called "ref- 
erence," system H' is constructed, which must be defined 
on the same lattice and must have the same interaction 
part as the original system H. Thus both the self-energy 
functional of the original system f2[Z] and the one of 
the reference system f2'[Z] contain the same F[Y], which 
can be eliminated by comparison from the expressions of 
the two self-energy functionals. This yields for bosonic 
systems^ 3 . 

0[Z] = n'[Z] - Tr ^(-(G'o -1 - Z)) 

+ Trln(-(G - 1 -Z)), (3) 

where quantities with prime correspond to the reference 
system and Go is the free Green's function. The free 
Green's function is defined as Gq 1 = (u + fi)! — T, where 
T contains the hopping matrix and all other one-particle 
parameters of the Hamiltonian except for the chemical 
potential /i, which is already treated separately in the 
definition. The symbol Tr denotes a summation over 
bosonic Matsubara frequencies and a trace over site in- 
dices. The self-energy functional 0[Z] given by Eq. ([3]) is 
exact. In order to be able to evaluate the functional, the 
search space of the self-energy Z has to be restricted^ 
which consists in an approximation. More precisely, the 
functional f2[Z] is evaluated for the subset of self-energies 
available to the reference system H' , see Fig.[T](a). Prac- 
tically, this is achieved by varying the single-particle pa- 
rameters of the reference Hamiltonian in order to find 
the stationary point of the grand potential. Thus the 
functional 0[Z] becomes a function of the set x of single- 
particle parameters of H' 

fi(x) = n'(x) - Tr ln(-(Go- 1 - Z(x))) 
+ Trln(-(G - 1 -Z(x))) 
= n'(x) + Tr ln(-G'(x)) - Tr m(-G(x)) (4) 

leading to the stationary condition 

^ = 0. (5) 

OX 

In VCA the reference system is given by the decom- 
position of the total system into identical clusters, see 
Fig.[T](b). In order to implement the Q-matrix approach, 
we solve each cluster by means of the band Lanczos 
method i 26 ' 35 The Green's function of the total system 
is obtained via the relation 

G- 1 (w) = G'- 1 (w)-V, (6) 
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(a) 



(b) 



FIG. 1: (Color online) (a) The search space for the self-energy 
T. is restricted to self-energies H(x) which are accessible via the 
reference system H' . (b) Lattice decomposition of a square 
lattice into 2x2 site clusters. 



which can be deduced from the Dyson equation of the 
total system G _1 = Gq 1 — Z(x) and the reference system 
G /_1 = Gq -1 — E(x). The self-energy can be eliminated 
and it follows that 

G- 1 = G'- 1 -(G^ 1 -G - 1 ). 

The expression in parenthesis defines the matrix 

V = G^ 1 - G^ 1 = ((w + - T') - ((w + /z)i - T) 

= -0i- M ')i + (T-T'). (7) 

With Eq. ^ the grand potential f2(x) can be rewritten 



n(x) = fl'(x) + Trln(i - VG') 



(8) 



The decomposition of the TV-site lattice into clusters of 
L sites can be described by a superlattice. The original 
lattice is obtained by attaching a cluster to each site of 
the superlattice.— A partial Fourier transform from su- 
perlattice indices to wave vectors k, which belong to the 
first Brillouin zone of the superlattice, yields the total 
Green's function 



G-^k, w) = G'- 1 ^) - V(k) 



(9) 



Due to the diagonality of G' in the superlattice indices 
its partial Fourier transform does not depend on k. The 
matrices in Eq. ([9]) are now defined in the space of cluster- 
site indices and are thus of size L x L. The TV wave 
vectors k from the Brillouin zone of the total lattice can 
be expressed as 



k + K 



(10) 



where K belongs to both the reciprocal superlattice and 
the first Brillouin zone of the total latticed 



B. Q-matrix formalism for bosonic systems 

The frequency integration implicit in the expression for 
the grand potential, given in Eq. ©, can be carried out 
analytically, yielding at zero temperatur e) 32 ' 33 ' 36 

n(x) = Q'(x)+J2K-^Yl E Mk), (ii) 



A' <0 



where A^. and A r (k) are the poles of the cluster Green's 
function and total Green's function, respectively. The 
number of clusters N/L is denoted as N c . The poles X' r 
of the cluster Green's function can be readily obtained 
from the Lanczos method, whereas the poles of the to- 
tal Green's function A r (k) can be evaluated with the so- 
called Q-matrix formalism, which was originally proposed 
for fermionic Green's functions! 26 ' 27 Here, we extend this 
formalism to the generic case, i.e., we include bosonic 
Green's functions. As we will see, this extension is non- 
trivial, since it involves non-unitary transformations. 

For zero temperature, the cluster Green's function 
reads^I 



-e(V>o|<4 



1 



w-(#'-wo) J 

1 

J oj + (H> - Lj Q ) 



a] |Vo> 
Oi\1>o) , (12) 



where |-0o) is the ground state of the N p particle system, 
wo is its (grand-canonical) energy, and e = 1 (e = — 1) for 
bosonic (fermionic) Green's functions. The first term on 
the right-hand side of Eq. (|12[) describes single-particle 
excitations from the N p particle ground state and can 
thus be referred to as particle term, whereas the second 
part corresponds to single-hole excitations and can be 
called hole term. Inserting the identity 1 = J2-y It) (tI 
into each part of Eq. (fl^)) . where I7) are the eigenvec- 
tors of the reference Hamiltonian with corresponding 
eigenvalues , yields the Lehmann representation of the 
Green's function 



G'ii M = E 



(ipo\cn I a) (a| a] |Vo) 
u-(u/ a - oj a ) 

g OJ + [U'f, - OJQ 



which can be cast into the form 

G'ii ( w ) = E ^'T — ^TT'^VyQ 



(13) 



(14) 



In Eq. (fl"4"|) , we have introduced the following notation 
nt - / (7l a) |Vo) |7> G H Np+1 



10 It) \i)£K Np 



and 



a: 



u\ - uj |t) S W-n p +i 



N v -1 



w - w 7 



l7>e« 



d^> 17) g n Np +i 



|7> e U Np -i 



(15) 



(16) 



(17) 



where T-Lm is the Hilbert space of an M particle system. 
With 



k A r (k)<0 



3 77' 



(18) 
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the cluster Green's function can be written in matrix 
notation 



G' = Qg' H SQt 



(19) 



With the help of this expression the VCA Green's func- 
tion Eq. © can be rewritten as 

G ^ G, T-W = Qg,SQf l-VQVsQt 
= Qg'SQ t {l + VQg'SQ t + ...} 



Qg'[l-SQtVQg'] 'SQt 
Q g'- 1 -SQtVQ SQT ' 



(20) 



where in the third step we expanded the fraction in a 
Taylor series. The matrix g' is diagonal and contains the 
poles of the cluster Green's function G', see Eq. (fT5|) . It 



can be written as g' 1 
Plugging this into Eq. ([20|) yields 

G = Q- 



A with A 



77' 



7 ^77 



-SQ 



(21) 



u- (A + SQtVQ)' 

We introduce the matrix M = A + S V Q. This matrix 
can be diagonalized as MX = XD, where D is a diagonal 
matrix containing the eigenvalues of M and X is the ma- 
trix of the eigenvectors of M. The eigenvalue equation 
of the matrix M can be rewritten as M = XDX -1 , where 
X" 1 is the inverse of X and not its transpose as M is a 
non-symmetric matrix. From that we obtain 

(uj-M)- 1 = X(uj -D)-^" 1 



Therefore, the poles of the total Green's function G in 
Eq. (|21[) are the eigenvalues of the matrix M . The ma- 
trices G and V are defined on the space of cluster-site 
indices. Thus G and V are of size L x L and depend on 
the wave vector k, see Eq. ©. The matrix Q is of size 
L x K, where K is the dimension of the Krylov space 
generated in the band Lanczos method. Due to the de- 
pendence of V on k the diagonalization of the matrix M 
yields K eigenvalues D rr i = X r (k)S rr i, which are used 
in Eq. (|lll) . The diagonalization has to be repeated for 
all wave vectors k. With that the grand potential f2(x) 
can be evaluated. The crucial point is that for bosonic 
Green's functions, the entries of the diagonal matrix S 
can be both 1 as well as —1, see Eq. (|17p . Therefore, the 
eigenvalue problem is not symmetric^ 

The factorization of the total lattice into clusters 
breaks the translational symmetry of the lattice. Hence 
the total Green's function would depend on two wave 
vectors k and k', which is certainly not correct for a 
periodic lattice. This has to be circumvented by a peri- 
odization prescription that provides a total Green's func- 
tion G(k, lu) depending only on one wave vector k. The 
periodization prescription proposed in Ref. [29| (Green's- 
function periodization) reads as follows: 



G(k, ^) = ^E e " 



* k ^-^G^(k, 



(23) 



where k is a wave vector of the total lattice and r Q refers 
to lattice sites a of the cluster. The wave vectors k in 
Eq. (|23|) can be replaced by the total wave vectors k as 
they just differ by a reciprocal superlattice wave vec- 
tor, see Eq. ([10]) ■ With Eqs. ([21]) and ([22]) the periodized 
Green's function can be rewritten in matrix notation 

G(k, u) =v[QX(w- D)- 1 X" 1 SQ t v k , (24) 

where the vector v k and its adjoint contain L plane 
waves 



1 



(e~ 4kro , e~ 4kri 



There exists as well an alternative periodization pre- 
scription where the self-energy Y. is periodized^ This 
self-energy periodization should prevent spurious gaps, 
which arise in the spectral function. However, at least for 
fcrmion systems, this procedure yields spurious metallic 
bands in the Mott phase for arbitrarily large U. Since we 
do not observe any spurious gaps in the spectral function 
of the 2D BH model we use the periodization on the 
Green's function defined in Eq. (1231 . 

With the wave- vector resolved Green's function of the 
total system G(k, w) we are able to calculate the single- 
particle spectral function 



A(k, u) = --ImG(k, u) , 

7T 



(22) the density of states 



N(w) 



f A(k, u,)dk=±Y, A Q*> 



(25) 



(26) 



and the momentum distribution 

n(k) = - ( A(k, w) duj . 

J — oo 

The frequency integration can be evaluated directly by 
means of the Q-matrix formalism, which yields a sum of 
the residues of the Green's function, see Eq. ([2~4"]) . corre- 
sponding to negative poles A r (k) < 0, 

»(k)= Y, (v k QX) r (X- 1 SQ t v k ) r . (27) 

A r (k)<0 

IV. RESULTS 

The BH model exhibits a quantum phase transition 
from a Mott to a superfluid phase when the ratio between 
the hopping strength and the on-site repulsion t/U is in- 
creased or when particles are added to or removed from 
the system. The Mott phase is characterized by an in- 
teger particle density, a gap in the spectral function and 
zero compressibility^ 



5 



2- 
1.5 

L 1, < 

0.5 

V 





□ L=2x1 




* L=2x2 


. * 


L=8 


* * e „ 




0.02 0.04 


0.06 


t/U 




(a) 






0.02 0.04 0.06 
t/U 

(b) 



FIG. 2: (Color online) Phase boundaries of the Mott phase 
of the 2D BH model (Mott lobes), (a) Results of our VCA 
calculation with various cluster sizes for the reference sys- 
tem. The geometry of the 8-site cluster is visualized in the 
inset. The gray shaded area indicates the results of the pro- 
cess chain approach . 21 ' 22 (b) Phase boundaries obtained for 
the 8 site cluster. The marks refer to the parameters where 
we evaluated the spectral functions. 



The first two Mott lobes of the 2D BH model obtained 
by means of VCA are shown in Fig.[2j We used the chem- 
ical potential x = {fj,} as variational parameter, which en- 
sures a correct particle density of the total system . 33 ' 39 
In contrast to the one-dimensional result a 33 ' 40 the Mott 
lobes of the 2D BH model are round shaped. The 
gray shaded area in Fig.[2](a) presents the phase bound- 
aries calculated within the process chain approach by 
N. Teichmann et al. in Refs. [U and H2, which arc ba- 
sically identical to the QMC results by B. Capogrosso- 
Sansone et al., see Ref. 1 14. The agreement is quite good 
for small hopping. However, VCA seems to overestimate 
the critical value of the hopping {t/U) c , which determines 
the tip of the Mott lobe. For the critical hopping of the 
first Mott lobe, we obtain approximately (t/U)]. — 0.067 
and for the second one (t/U) 2 = 0.038. Latest process 
chain approach j 21 ' 22 QMC (Ref. Il4| ) and strong-coupling 
perturbation theory^ 9 - results yield (t/U)] = 0.059 and 
(t/U) 2 — 0.035 for the critical parameter of the first and 
second Mott lobe, respectively. 

The spectral functions A(k, uj) and the densities of 
states N(u>) for parameters of the first Mott lobe are 
shown in Fig.[3J The spectral function is displayed on the 
conventional path around the Brillouin zone k = (0, 0) 
over (it, 7r) to (tt, 0) and back to (0, 0), and we use an ar- 
tificial imaginary-frequency broadening 77 = 0.05. A pe- 
culiarity of bosonic systems is that the hole band of the 
spectral function has negative spectral weight whereas 
the particle band has positive spectral weight. This fol- 
lows from the definition of the bosonic Green's function 
which has a negative sign in front of the hole term, sec 
Eq. (|12[) . In the figures we always plot the absolute value 
of the spectral function. The local density of states 
is defined as a wave- vector summation of A(k, ui), see 
Eq. (|26p . Therefore we observe a negative peak in the 
density of states, which corresponds to the hole band 
of the spectral function. For bosonic Green's functions 
the density of states is not a probability distribution, as 
it contains negative values. Taking the absolute value 
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FIG. 3: (Color online) Spectral function A(k, uj), left column, 
and density of states N(lj), right column, in the first Mott 
lobe for the parameters (a) t/U = 0.005, fi/U = 0.5, (b) 
t/U = 0.03, n/U = 0.4 and (c) t/U = 0.06, fi/U = 0.35. The 
captions of the subfigures refer to the marks in Fig.[5](b). 



would yield an all positive density of states, however, it 
would not be normed and is thus no probability distri- 
bution either. For small hopping, the gap in the spectral 
function is large and the bands are rather flat, i.e., the 
width of the bands is small, see Fig. [3] The correspond- 
ing density of states contains two well-separated peaks. 
For increasing hopping, the gap of the spectral function is 
decreasing and the width of the bands is increasing. Pur- 
suant to the spectral function, the peaks in the density 
of states become broader for increasing hopping. The 
intensity of the two bands is almost constant for small 
hopping independent of the wave vector k, whereas for 
large hopping a large intensity can be observed at k = 0. 

The boundaries of the Mott lobes correspond to the 
chemical potential of the state with one additional par- 
ticle (hole), which is obtained directly from the single- 
particle (single-hole) minimum excitation energy. For 
this reason, we evaluate the phase diagram in Fig. [2] by 
taking the minimal gap of the spectral function for each 
t/U , which always occurs at k = 0. 

The spectral functions and densities of states in the 
second Mott lobe corresponding to the marks IV and V 
in Fig.[5](b) are shown in Fig. 2] Qualitatively they are 
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FIG. 6: (Color online) Direct comparison (a) between the 
momentum distribution n(k) obtained by means of VCA and 
a strong-coupling perturbation theory with scaling ansatz.— 
The data identified with the letters SC correspond to the 
strong-coupling results. The relative deviations between VCA 
results and strong-coupling results with scaling ansatz are 
shown in (b). The Roman numerals in the legends refer to 
the parameters marked in Fig.[2](b). 



FIG. 4: (Color online) Spectral function A(k, u), left column, 
and density of states N(lo), right column, in the second Mott 
lobe for the parameters (a) t/U = 0.01, fi/U = 1.5 and (b) 
t/U — 0.035, fi/U = 1.4. The captions of the subfigures refer 
to the marks in Fig.[2](b). 
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FIG. 5: (Color online) Momentum distribution n(k) in (a) the 
first Mott lobe and (b) the second Mott lobe. The Roman 
numerals in the legends refer to the parameters marked in 
Fig-i(b). 



very similar to the spectral functions and densities of 
states in the first Mott lobe. Particularly, the intensity 
distribution of the bands seem to be strongly related. Yet 
the peaks of the density of states are larger due to the 
twice as large particle density within the second Mott 
lobe and thus the absolute value of the spectral weight 
in the second Mott lobe is larger than the one in the first 
Mott lobe. 

The momentum distribution n(k) corresponding to the 
spectral functions in the first and second Mott lobe are 
shown in Fig.[5j The particle density in the first Mott 
lobe is one and thus n(k) is centered around one in 
Fig.[5](a). For the second Mott lobe n(k) is centered 
around two, see Fig.[5](b). The particle density n(k) is 
extremely flat for small hopping whereas it is peaked at 
k = for large hopping, which is already a precursor for 
the Bosc-Einstcin condensation where all particles con- 
dense in the k = state. This behavior directly reflects 



the intensity distribution of the bands in the spectral 
function. There is excellent quantitative agreement be- 
tween our VCA results for the momentum distribution 
and results obtained by means of QMC and a strong- 
coupling perturbation theory with scaling ansatz,— see 
Fig. [6] We compare the momentum distributions for the 
parameters I (t/U = 0.005) and II (t/U = 0.03), and 
observe that the relative deviations between our VCA 
results and an approach obtained by combining strong- 
coupling perturbation theory with a scaling ansata^ are 
almost zero for small hopping t/U = 0.005 and less than 
one percent for medium hopping t/U = 0.03. This lat- 
ter methods is certainly more accurate than VCA in the 
evaluation of the momentum distribution. However, it 
should be mentioned that the information about the criti- 
cal point (critical exponents and critical hopping strength 
(t/U) c ) have to be inserted "by hand," in order to op- 
timize the results. This information, in turn, must be 
extracted, e. g., from a QMC calculation. On the other 
hand, our VCA results are obtained directly without the 
need to introduce external parameters. 



V. CONCLUSIONS 

In the present paper, we presented and discussed re- 
sults obtained within the variational cluster approach 
for the spectral properties of the two-dimensional Bose- 
Hubbard Hamiltonian. This is a minimal model to 
describe bosonic ultracold atoms confined in optical 
lattices^ and it undergoes a quantum phase transition 
from a Mott to a superfluid phase depending on the 
chemical potential (i, and the ratio between the hopping 
strength and the on-site repulsion t/U. In particular, 
we determined the first two Mott lobes of the phase di- 
agram and found reasonable agreement with essentially 
exact results from QMC simulations and from the process 
chain approach. In particular, the variational cluster ap- 
proach yields very good results for the phase boundaries 
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apart from the region close to the lobe tip. Here, strong- 
coupling expansions and QMC calculations are, clearly, 
much more accurate. Yet it should be emphasized that 
the computational effort is considerably lower for VCA 
than for QMC. Furthermore, we evaluated spectral func- 
tions in the first and second Mott lobe. An important 
aspect of VCA is that the Green's function of the sys- 
tem is obtained directly in the real frequency domain, 
which allows for a direct calculation of the spectral func- 
tion. On the other hand, QMC quite generally provides 
correlation functions in imaginary time. Imaginary-time 
correlation functions have to be analytically continued to 
real frequencies, which is a very ill-conditioned problem, 
as the data contain statistical errors. In QMC this ana- 
lytical continuation is best carried out by means of the 
maximum entropy method. A very accurate dispersion 
(without spectral weight) has been also obtained by a 
strong-coupling expansion^ The intensity distribution 
of the spectral weight is similar for the spectral functions 
of both Mott lobes, leading to an evenly distributed spec- 
tral weight for small hopping strengths and to a distribu- 
tion sharply peaked at k = for large hopping strengths. 
The latter indicates a precursor to the Bosc- Einstein con- 
densation occurring above a certain critical hopping. We 
also evaluated the densities of states and momentum dis- 
tributions corresponding to the calculated spectral func- 



tions. We compared our VCA results for the momentum 
distribution with strong-coupling perturbation-theory re- 
sults, where a scaling ansatz has been used, and found 
excellent quantitative agreement. Finally, as a techni- 
cal point, we extended the Q-matrix formalism to deal 
with bosonic Green's functions, which, in contrast to 
the fermionic case, produces a non-symmetric eigenvalue 
problem. 
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